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We consider a voltage-biased nanoelectromechanical Josephson junction, where a suspended 
nanowire forms a superconducting weak-link, in an inhomogeneous magnetic field. We show that a 
nonlinear coupling between the Josephson current and the magnetic field generates a Laplace force 
that induces a whirling motion of the nanowire. By performing an analytical and a numerical anal- 
ysis, we demonstrate that at resonance, the amplitude-phase dynamics of the whirling movement 
present different regimes depending on the degree of inhomogeneity of the magnetic field: time 
independent, periodic and chaotic. Transitions between these regimes are also discussed. 
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I. INTRODUCTION 

Over the past fifteen years, nanoelectromechanical sys- 
tems (NEMS) have been widely used for the exploration 
of the quantum world and for the development of new 
technological applications [1-7]. Because of their size, 
high frequency operation, small mass and high perfor- 
mance, NEMS are currently considered promising candi- 
dates for achieving the quantum limit of mechanical mo- 
tion. It is expected that the quest for the quantum regime 
in such devices will elucidate questions of fundamental 
nature in physics, for instance, the quantum-mechanical 
description of macroscopic objects [8-13]. Similarly, 
physical systems involving NEMS resonators are excel- 
lent tools for theoretical and experimental studies of 
nonlinear dynamical systems [14, 15]. Indeed, examples 
of complex dynamical phenomena in NEMS are numer- 
ous and include chaotic behavior [16-19], bifurcation- 
topology amplification [20], nonlinear switching dynam- 
ics [21] and nonlinear frequency pulling [22] to name but 
a few. 

By making use of the potential offered by NEMS res- 
onators, G. Sonne et al. [23] studied the nonlinear dy- 
namics of a suspended carbon nanotube coupled to two 
voltage biased superconducting electrodes. In presenting 
their work, Sonne and collaborators assumed that the 
nanoelectromechanical junction was subjected to a ho- 
mogeneous magnetic field perpedicular to the axes of the 
nanowire. For such a system, the authors demonstrated 
the possibility to pump energy from the electronic sub- 
system into the mechanical vibrations; they also demon- 
strated that the system had more than one regime of 
finite-amplitude stationary nonlinear oscillations. In par- 
ticular, a region of bistability was found and the authors 
showed that it should be detected in the corresponding 
dc Josephson current (see discussion in Rcf. 23). 

In this article, we consider the same voltage-biased 
nanoelectromechanical system studied by Sonne and co- 
workers [23], but now we extend the analysis to the case 
in which the NEMS resonator is subjected to a nonuni- 
form magnetic field. As will be discussed below, inhomo- 



geneity of the field causes the conducting nanoresonator 
to execute a whirling movement resembling a jump-rope 
like motion. The purpose of this paper is to analyze 
the time evolution of the amplitude and relative phase 
of the nanotube whirling motion. We will demonstrate 
that the coupled amplitude-phase dynamics exhibit dif- 
ferent stationary regimes depending on the degree of the 
magnetic field inhomogeneity: time independent, peri- 
odic and chaotic. 

The scope of this article is parallel to research con- 
ducted by Conley et al. [24] and Chen et al. [25], who pre- 
viously investigated nonlinear and nonplanar dynamics 
of suspended nanowires excited by an electrostatic force. 
Here, however, we consider another type of driving mech- 
anism: the Laplace force. This excitation mechanism to 
accomplish whirling motion of suspended nanoresonators 
has not been considered before. 

The article is organized in the following manner. In 
Sec. II we introduce the model Hamiltonian for the con- 
sidered system and examine the coupling between the 
mechanical and electronic subsystems. At the end of the 
section, we derive the equations of motion governing the 
dynamics of the nanoresonator driven by an ac Josephson 
current. In Sec. Ill we obtain the differential equations 
for the time evolution of the amplitudes of the nanotube 
vibration and discuss their dynamics through the analysis 
of the numerically computed results. Finally, in Sec. IV 
we provide the summary and concluding remarks of this 
work. 



II. MODEL HAMILTONIAN SYSTEM 

The diagram in Fig. 1 is a schematic illustration of a su- 
perconducting hybrid nanostructure, a superconducting- 
normal-superconducting (S-N-S) nanoelectromechanical 
Josephson junction driven by a dc voltage bias V. The 
junction consists of a metallic carbon nanotube of length 
L suspended between two voltage-biased superconduct- 
ing leads. In such a geometry, the nanotube is simulta- 
neously serving as a mechanical resonator and as a weak 
link between the superconducting electrodes. In the lay- 
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FIG. 1: (Color online) (a) Schematic diagram of the voltage- 
biased S-N-S nanoelectromechanical Josephson junction con- 
sidered in the article. A doubly clamped metallic carbon nan- 
otube suspended over a trench of length L, forms a weak link 
between two voltage-biased superconducting electrodes. The 
junction is influenced by an inhomogeneous magnetic field H, 
generated by a wedge-shaped MFM cantilever tip at a dis- 
tance d from the nanotube at rest (blue circle), (b) Nanotube 
displacement (green circles) in the z-y plane induced by a 
Laplace force. Rt is the curvature radius of the magnetic tip. 

out of the system, the NEMS junction is under the influ- 
ence of an external inhomogeneous magnetic field H, gen- 
erated by a magnetic force microscope (MFM) cantilever 
tip parallel to the axis of the nanotube at a distance d. 

As demonstrated by Sonne et al. [23], the electronic 
subsystem can be effectively modeled in terms of the 
Bogoliubov-de Gennes Hamiltonian 

n = J ^\x)(n + n A )^(x)dx, (i a ) 

h 2 f d . ieA x \ 2 A , . . , . 

Wo = -^(^-^-rJ + ° Mxh (lb) 

Ha = A(x) [a x cos <j>(t) + sgn(x)a y sin <p(t)] . (Ic) 

Here, 4/(x) and its complex conjugate 4* (x) are two- 
component Nambu spinors that annihilate and cre- 
ate electrons/holes at position x, respectively. In the 
Bogoliubov-de Gennes equation, Ho is the Hamilto- 
nian for electrons/holes in the metallic carbon nanotube, 
where U(x) is the potential barrier between the tube 



and the bulk superconducting leads, A x is the mag- 
netic vector potential in the x-direction and <Xj are the 
Pauli matrices. The superconducting electrodes are de- 
scribed by Ha, where the gap function is given by 
A(x) = AoO(2|x| — L), with Ao the order parameter and 
0(x) the Heaviside step function. The phase difference 
between the superconducting leads is denoted by cj> and, 
in accordance to the second Josephson relation, its time 
evolution is <p(i) = u>jt, with wj = 2eV/h the Josephson 
frequency. 

We now consider a magnetic field of the form H = 
(0, H y , H z ), which can be generated, for instance, by an 
MFM cantilever tip in the form of a wedge, cf. Fig. 1. A 
first order Taylor expansion of the magnetic field yields 
H = {0,yd y H y (x,0,0),H z (x,0,0)+zd z H z (x,0,0)), 
where H z (x,0, 0) and diHi(x, 0, 0) represent the magni- 
tudes of, respectively, the magnetic field z-component 
and the magnetic field gradients, both evaluated at 
the axis of the nanotube. A straightforward calcula- 
tion from the Maxwell equation V • H = 0, indicates 
that d y H y (x,0,0) = -d z H z (x,0,0) = -H' z . Setting 
H z (x, 0, 0) = H z , the magnetic field reads 

U = (0,-H' z y,H z +H' z z) . (2) 

This field is also obtained from the equation H = V x A, 
where the vector potential A is given by 

A = -{[H z + H' 8 z] y, 0,0). (3) 

By substituting the magnetic vector potential, Eq. (3), 
into the Hamiltonian given by Eq. (Ib), we can read- 
ily observe that Ho is altered by deflecting the nanotube 
in the z-y plane. In this case, non-planar whirling dis- 
placement of the nanotube is due to the Laplace force 
generated by coupling between the magnetic field and 
the Josephson current. As a consequence of the motion 
in the magnetic field, an electromotive force is induced 
along the nanomechanical weak-link and its magnitude 
depends on the rate of change of the nanowire profile in 
the z-y plane, i.e., the rate of change of the magnetic flux 
through the circuit. In this description, now, the super- 
conducting phase difference is not only a function of the 
bias voltage, but also of the nanowire deflection in the z-y 
plane. We decompose the nanotube motion in this plane 
into two independent deflections y(x,t) — uo(x)a(t) and 
z(x,t) = uo(x)b(t), where Uq(x) is the normalized and 
dimensionless profile of the fundamental mode in both 
directions. Then, following Shekhter et al. [26], a unitary 
transformation which moves the vector potential from 
the kinetic part of Ho to the phase difference between 
the superconducting electrodes is applied. As a result, 
the expression for the superconducting phase difference 
has the form 

^ t) = m + ^L [hz+Km]a{th (4) 

The parameters h z — aH z and h' z — j3H' z are the renor- 
malized magnetic field and magnetic field gradient in the 
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FIG. 2: (Color online) Numerical simulations of Eqs. (6a) and (6b) for the time averaged coordinates Y and Z as function of 
the driving voltage V. The plots are calculated for 7 = 0.001, e = 0.012, h = 1 (L = 1/im, h z = 40 mT, Rt ~ 66 nm). 



z-direction calculated at the axis of the tube, respectively, 
with a, /? ~ 1 correctional factors originating from geo- 
metrical considerations. 

The nanotube mechanics is thus described through the 
projection amplitudes a(t) and b(t). For the conjugate 
variables {a(t),p a (t)} and {b(t),pb(t)} (j)j(t) denoting 
the generalized momenta) one can formulate the follow- 
ing Hamiltonian function 

H{p a ,p b , a, M) = ^ (p2 + Pb) + ^y- + b 2 ) - 

- 2DA co$(ip(a,b, t)) , (5) 

where m and ui are the mass and the mechanical eigen- 
frequency of the nanoresonator, respectively. The last 
term in Eq. (5) corresponds to the Josephson energy 
Ej(if(a,b,t)), with D the transmission coefficient of the 
junction. The equations of motion for a(t) and b(t) 
are then obtained from the Hamilton equations. Writ- 
ten in terms of the dimensionless deflection coordinates 
Y(t) = [AeLh z /h]a(t) and Z(t) = [AeLh z /h]b(t), the 
resulting set of differential equations for the nanotube 
amplitudes is: 

Y(t)+jY{t) + Y(t) = 

- £ [1 + *Z(t)] sin (Vt + [1 + hZ(t)] F(t)) , (6a) 

Z(t)+jZ{t) + Z(t) = 

- exY{r) sin (ft + [1 + kZ(t)\ Y(r)j . (6b) 

Here, we have added a dimensionless phenomenological 
damping coefficient 7 = [7/mw]. In these equations 
e = [8eL 2 hlj c /mHu} 2 ] with j c = [DA e/2h], the critical 
current through the junction. We also set the timescale 
to t = u)t and, consequently, V = [eV/fku] . The parame- 
ter 

H = AeLh z R T ' ^ 



where Rt = \h z lh! z ^ denotes the curvature radius of the 
magnetic cantilever tip, characterizes the degree of inho- 
mogeneity of the magnetic field and will be referred to 
as the control parameter. It turns out that x determines 
the dynamical behavior of the system. One can then re- 
alize the significance of x by considering fixed system 
parameters (L,h z ) and by letting the Rt vary. In the 
limit Rt — > 00, the control parameter vanishes and the 
equations of motion given by Eqs. (6) reduce to the case 
discused by Sonne et al. [23] where the magnetic field 
is uniform and, therefore, the amplitude of the driving 
force acting in the z-direction becomes zero. 



III. NUMERICAL RESULTS AND DISCUSSION 

For a qualitative and quantitative discussion of the dy- 
namic behavior of the nanowire amplitudes, we consider 
the following system parameters: a carbon nanotube of 
radius r = 1 nm and length L = 1 fj,m [27], supercon- 
ducting order parameter Ao ~ 10 meV, Josephson critical 
current j c ~ 100 nA [28], and quality factor Q — 10 3 [29]. 
We also assume that 7 = 1/Q, h z — 40 mT, e = 0.012, 
and Rt is varied from 53 nm to 544 nm. 

Numerical simulations of Eqs. (6a) and (6b) allow us to 
study the time average of the nanotube deflection coor- 
dinates Y(t) and Z(t) as functions of the driving voltage 
V . In doing so, one can notice that the system response 
presents a series of resonance peaks at integer values of 
the driving voltage, i.e., the amplitude of the nanores- 
onator is not damped provided the resonant condition is 
fulfilled: the Josephson frequency matches the mechan- 
ical frequency (see Fig. 2). This resonant phenomenon 
was first studied by Sonne et al. [23], who attributed a 
direct resonance at V = 1 and a parametric resonance at 
V = 2. Accordingly, the same conclusion can be drawn 
from the results presented in Fig. 2. In the remainder of 
the article we will be mainly focusing on the parametric 
regime and take V = 2. In this case, the dynamic be- 
havior of the amplitudes can be analized by postulating 
a solution for both deflection coordinates in the form 
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FIG. 3: (Color online) Numerical simulations for the stability analysis of Eqs. (9). Stable and unstable stationary points are 
indicated by colored diamonds (0) and black crosses ( x ) respectively. Diamonds of identical color indicate the four envelopes 
of the stable solution X s , crosses in the same region (I or II) of both phase planes belong to the same unstable solution X u . 
The zero solution, X = 0, is unstable. Here, 7 = 0.001, e = 0.012, x = 0.12 (L = 1 fj,m, h z = 40 mT, R T = 544 nm). 



Y(t) = A r (r) cos(t) + Ai(r) sin(r) , (8a) 
Z(t) = B r (r) cos(r) + B^t) sin(r) . (8b) 

On condition that 7, e, ex <C 1, the four envelopes in the 
vector X = ( A r (r), Aj(r), B r (r), Bi(r)) vary slowly in 
time, i.e., dX/dr <C 1 and, an averaging method [30] can 
be employed in order to derive the equation of motion for 
X. By substituting the ansantz provided in Eqs. (8) into 
the system of equations in Eqs. (6) and integrating over 
the fast oscillations one gets 

dA r jA r dg dB r ^B r dg 

~dV + = dAi ' lh r + ~T~ = dBi ' . , 
dA t 7A: _ dg dB, jB, _ dg ^ ' 

~dV + ~Y ~ ~d~A^' "d7 + ^~~~aB^' 

where 

g(x) = / cos([l + hB t cos(9) + ycBi sin(9)] 

27r J-TT 

x [A r cos(&) + A iS hx(&)+2@])de. (10) 

Here, g{x) is the generatrix Hamiltonian function. The 
study proceeds by performing a stability analysis based 
on Eqs. (9). In general, we shall study solutions of a sys- 
tem of coupled ODEs, dX/c?r = f(X, x), and solutions 
of a system of algebraic equations = f(X, x). 

The discussion commences by highlighting the sym- 
metric nature of the stationary solutions for the algebraic 
system in Eqs. (9). For the above considered system pa- 
rameters with Rt = 544 nm and k = 0.12, numerical 
simulations for the slowly-varying envelopes show that 
the system presents several stationary points. Such qual- 
itative behavior is visualized in the phase space diagrams 
for the nanotube amplitudes, Y and Z, in Fig. 3. There, 



solutions pictured by colored diamonds (0) and black 
crosses ( x ) correspond to stable and unstable stationary 
points, respectively. Similarly, the phase portraits are di- 
vided in two regions, denoted I and II, by diagonal sym- 
metry axes. The symmetry axes in the Y and Z coordi- 
nates satisfy the relations Ai — A r = and Bi + B r = 0, 
respectively. From Fig. 3 it follows that the presented 
pattern of solutions has a mirror-image symmetry with 
respect to the line of symmetry that divides the phase 
diagrams into regions I and II, where diamonds of the 
same color represent the four components of a stable com- 
puted solution X s and crosses in the same region (I or 
II) of both phase spaces account for the envelopes of an 
unstable solution X M . From a symmetric point of view, 
the system is characterized by two stable and two unsta- 
ble (including X = 0) stationary points in each region. 
However, the dimensionless angular momentum of the 
nanowire, "L x = YdZ/dr — ZdY/dr, in both regions is 
equal in magnitude but opposite in sign. 

Continuing our exploration, the computed results in- 
dicate that the bistable regime is only attained for x < 
0.120. By letting the control parameter increase further, 
bistability is abandoned and the number of stationary so- 
lutions of Eqs. (9) is modified as well as their stability. In- 



TABLE I: Bifurcation pattern in the symmetric regions of 
solutions I and II as a function of the control parameter >c. 
Solution X = is included. 





Number of Number of 


x-interval 


stable solutions unstable solutions 


x < 0.120 


2 2 


0.120 <x < 0.125 


3 3 


0.125 <h< 0.140 


1 3 
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FIG. 4: Period-doubling bifurcations for different values of the control parameter. Plots are shown for the first symmetry region 
of both phase spaces. Simulations were obtained for 7 = 0.001, e = 0.012 (L = 1 pim, h z — 40 mT). 



deed, Table I presents the number of stable and unstable 
stationary solutions in the interval 0.120 < x < 0.140 in 
both regions for these equations. As can be seen from the 
table, the system displays a phenomenon called branch- 
ing or bifurcation, which is a distinctive fingerprint of 
nonlinear dynamical systems [31]. The exhibited branch- 
ing pattern formation can be described as follows. For 
x < 0.120, bistability is accompanied by two unstable 
stationary points, one of them corresponds to X = (see 
Fig. 3). Once the control parameter is slightly increased 
over the threshold value x = 0.120, the nonzero unstable 
stationary point becomes stable, and simultaneously, two 
new unstable solutions appear and the system has now 
three stable and three unstable (including the trivial so- 
lution) stationary points. Next, by varying the control 
parameter in the interval 0.120 < x < 0.125, the two 
original stable points and the two unstable ones move in 
phase space and approach each other. Finally, when the 
control parameter ranges fom 0.125 to 0.140, the bistable 
and the created unstable solutions coalesce into two un- 



stable solutions. In this stage, there is only one stable 
and three unstable stationary points. 

It turns out that for x > 0.140, the system leaves 
the regime of equilibria and enters into the one of pe- 
riodic solutions. In fact, at the critical control parameter 
x ~ 0.150, there is an exchange of stability from the 
unique stable equilibrium to a stable limit cycle. This 
transformation in phase space is performed through a 
Poincare-Andronov-Hopf bifurcation and the computed 
results show that the stable limit cycle grows in phase 
space for 0.150 < x < 1. Hereupon, further changes in 
the control parameter will be reflected in the periodic- 
ity of the limit cycle as it can clearly be seen in Fig. 4. 
In this figure, the nanotube dynamics undergoes suce- 
sive period-doubling cascade bifurcations in the ampli- 
tude modulation when varying x. Due to the symmetric 
character of the solutions, the results are only plotted in 
region I. 

Period-doubling bifurcations pave the way for chaotic 
dynamics [32] and for the control parameter x ~ 1.25, 
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FIG. 5: (Color online) Chaotic whirling motion of the suspended carbon nanotube. Red and blue dots represent the dynamical 
flowing of the system for initial conditions in the first and second region of symmetry, respectively. In this case, 7 = 0.001, 
e = 0.012 , x = 1.25 (L = 1/im, h z = 40 mT, R T = 53 nm). 



a particular event takes place. In Fig. 5 the phase por- 
traits for the nanowire amplitudes suggest that the ini- 
tially disconnected symmetric regions I and II are now 
connected in a very complicated way. Furthermore, the 
flow diagrams for initial conditions in region I (red dots) 
and region II (blue dots) indicate that the system dy- 
namically evolves in the two regions, i.e., the flow lines 
generated by Eqs. (9) are sensitive to the initial values 
and are dense in both regions, where two strange attrac- 
tors can be readily identified. Hence, the system is in the 
chaotic regime. 

IV. CONCLUSIONS 

We have considered a voltage-biased nanoelectrome- 
chanical Josephson junction, where a suspended nanowire 
is serving as a weak link between two superconducting 
electrodes, in an inhomogeneous magnetic field. For our 
case study, we have assumed that the magnetic field is 
generated by an MFM cantilever tip and the nanowire 
is in the form of a metallic carbon nanotube. In such 
a scenario, the inhomegeneity of the field in conjunction 
with the Josephson current flowing through the tube, give 
rise to a Laplace force that induces the nanoresonator 
to perform a whirling movement. We have studied the 
time evolution of the amplitude and relative phase of this 
non-planar whirling motion and demonstrated that at the 
parametric resonance, their coupled dynamics exhibit a 
rich dynamical behavior characterized by: multistability, 
limit cycles and chaos. These stationary regimes depend 



on the degree of the magnetic field inhomegeneity, which 
in the present case, is related to the curvature radius of 
the magnetic cantilever tip. 

The experimental implementation of the system con- 
sidered in the article is plausible in light of current state- 
of-the-art nanofabrication techniques. For instance, J.- 
D. Pillet and collaborators designed and constructed a 
superconducting hybrid nanostructure that comprises a 
carbon nanotube suspended between two superconduc- 
tors [33]. Concerning the MFM cantilever tip, the ex- 
perimental results reported by Victor N. Matveev and 
co-workers in Ref. [34] suggest that it is possible to fab- 
ricate cantilever tips coated by a magnetic film with a 
curvature radius in the range of 50 — 70 nm with max- 
imum magnetic fields in the range of 40 — 80 mT. Due 
to a growing interest in complex behavior in nanodevices 
[35], the nonlinear and nonplanar phenomena exhibited 
by the system studied here have potential applications 
in signal processing, chaotic encryptation and random 
number generation [24, 25]. 
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